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Abstract 

Growth and folding in one-layered model tissue sheets are studied in a stochastic, 
lattice-free single cell model which considers the discrete cellular structure of the tissue, 
and a coarse grained analytical approach. The polarity of the tissue sheet is considered 
by a bending term. Cell division gives rise to a locally increasing metric. An exponential 
and a power-law growth regime are identified. In both regimes folding occurs as soon 
as the bending contribution becomes too small to compensate the destabilizing effect 
of the cell proliferation. The potential biological relevance is discussed. 

Many growth models have been inspired by biological systems or can be applied to them, 
see e.g. [1-5]. Most models consider either surface or bulk properties of more or less compact 
assemblies of particles (here: cells) as e.g. in Eden-like models of tumor- Q], bacterial- 0] or 
yeast growth 0. We here present a novel microscopic growth model where cell division can 
take place only within a single cell layer such, that a one-layered structure is maintained. 
In computer simulations with this model and in a coarse grained continuum approach it is 
shown that the growing layer cannot remain smooth but must fold by a mechanism that is 
believed to be generic, i.e. inherent to this structure. 

Such situations can occur in all biological cell systems in which cell division has to maintain 
a single- layered structure, e.g. in the skin 0,0] or in crypts, the cell proliferating units 
in the intestine. Strong perturbations of the physiological cell division in crypts result in 
strong folding producing patterns analogous to those observed during the formation 

of polyps or adenoma which constitute pre-patterns of intestinal cancer 0. 
For simplicity and in order to obtain a clear illustration of the underlying folding principle 
we focus on simple one-dimensional (lei) tissue manifolds ("cell chains") in two-dimensional 
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FIG. 1: (a) Cell division algorithm: During growth a cell deforms from a spherical shape into a 
dumb-bell by elongation of the axis from a = ^ a = 2R by small steps, (b) r2 is the local radius 
of curvature for cell i = 2, P2 its curvature angle. During deformation cell 2 attempts to orientate 
its axis into the direction of the local tangent to the circle in point P2- 

{2d) space. Our basic model assumptions are the existence of attractive nearest-neighbor 
(NN) interactions between cells to maintain the integrity of the (one-layered) tissue sheet, a 
bending energy that models the polarity of the cells in a sheet and cell division that takes 
into account potential size changes of the sheet. In order to maintain a one-cell thick struc- 
ture, cells have to be placed in the sheet again after cell division. Since cells grow and divide 
only if this does not result in too strong cell deformations or compressions, cell division 
at constant rate requires the migration of either the cell layer itself or of cells within the 
cell layer. The growth law and the geometry of the tissue sheet is then determined by the 
competition between the area proliferating and destabilizing cell growth and the stabiliz- 
ing bending energy of the tissue that locally confines cell movements perpendicular to the 
layer and thus can be responsible for a confinement of cell division. As soon as the bending 
energy becomes too small to smoothen local undulations that are stochastically created by 
local fluctuations in the proliferation of length, the layer roughens. If this occurs before cell 
deformations or compressions become strong enough to hinder cell division, the cells divide 
exponentially fast. Otherwise first the layer stays smooth but the growth law changes to 
sub-exponential growth. Later again, the layer folds. A large growth rate results in early 
crossovers. 

Below we first introduce the microscopic model that shows the crossover between both growth 
regimes as well as the folding. The instability within each growth regime is explained in an 
analytical approach afterwards. The basic unit in the microscopic model is a single cell. Each 
cell is assumed to be spherical directly after cell division and deforms during mitosis into a 
dumb-bell (Fig. la, |TD[), i.e. actively changes its 'equilibrium shape'. A particular config- 
uration of cells is characterized by its total energy, V^*°* = J2i<j ^ + J2i V^^^'^ + Y^i Vf^K 
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The terms on the rhs. characterizes (i) the NN-interaction, (ii) the bending contribution, 
and (iii) non-spheric cells. We assume that the cell-cell-NN, short-range potential has one 
minimum which results from the competition between (a) attractive interactions due to 
adhesion molecules anchored in the cell membranes and (b) repulsive contributions from 
the limited cell deformability and compressibility and the loss of membrane steric entropy. 
The NN-potential used in the simulations is = e {{{{2dij{t) - AR) /5) - if - 1) for 

2R < dij{t) <2R + 5 and V^^ = oo for dij{t) < 2R or dij{t) > 2R + 5 Here, dij{t) 
denotes the distance between the nearest circles of the neighboring dumb-bells i and j (a 
circular cell is a dumb-bell with axis length a = 0). The lower cutoff at dij = 2R takes 
into account the limited cell deformability and corresponds to an excluded volume. 6 de- 
termines the range over which a cell can be stretched or compressed in a certain direction. 
A typical value is 5 ~ 0.2R. e (> 0) is a measure for the resistance against deformations 
(e = 10 throughout this article). We introduce polarity by assuming that an anisotropic 
distribution of cell-cell adhesion molecules gives rise to a bending energy in addition to the 
isotropic energy contribution: V^"^"'''' = (1/3) Z]jii-i(^/2)((l/Ti) — Co)^rj/5j. n is the bending 
rigidity, cq the spontaneous curvature, the local radius of curvature and Pi the local angle 
of curvature (Fig. lb). Cell division in many tissue sheets is directed, e.g. during early 
embryogenesis, in the crypt and the skin, suggesting that cells are able to detect the position 



of their neighbors in order to determine the direction of their division . Accordingly we 
assume that non-spherical cells orientate their axis in a way that guarantees the cell division 
takes place in the direction of the tangent to the local radius of curvature (Fig. lb). This 
insures the maintenance of a one-cell thick tissue sheet. Formally this is carried out by an 
energy contribution Vl'°* = 7(0, — a°^*)^ that for 7^1 sufficiently penalizes deviations from 
the optimal orientation of the cell axis. Here, ctj is the momentary, the optimal angle 
of the axis in an external coordinate system. 

Active cell deformations during mitosis cause a pressure on the neighbor cells in the direc- 
tion of the deformation. This leads to an increase of the total energy \/*°*. We assume 
that the neighbor cells either move their center of mass or change their orientation in or- 
der to minimize V^*°*. Under the assumptions that (i) inertial terms are small compared to 
dissipative terms, (ii) processes not explicitly considered, such as the cell metabolism, intra- 
cellular movements of the cytoplasm and the reorganization of the cytoskeleton give rise to 
a stochastic component in the displacement of the cells, the dynamics can be modeled by 
the Metropolis method which corresponds to the numerical integration of a Master equation 
p|, |13|: Cells are randomly chosen to perform either a small translation step, orientation 
change or growth trial. Each translation or rotation is accepted with probability Pq = 1 if 
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FIG. 2: (a) Instability in a stretched cell configuration with fixed ends (snapshots, single arrows 
indicate the time direction), (a.l) start configuration, N = 100, k = 100 (a. 2): = 200, N = 110, 
(a. 3): = 50, N = 110. (b) Instability in a circular growing cell arrangement, time series for 
r-i = 1.25, K = 2000 for N = 8,16,32,64,128,256 (cf. Fig.3b). (c): snapshots for N = 512, 
r-i = 500, (c.l): k = 20, (c.2): k = 100 (cf. Fig.3a, points A and B). 

^ytot ^ ytot^^ _ ytot ^ q g^^^ ^j^j^ probability Pa = exp(-Ar*°VFT) if AV*"* > (hence 
isolated cells move diffusively in accordance with ref. [Q)- is the time between two 
migration or rotation trials (= 1 in the simulations). Ft is the fc^T-analog in cellular sys- 
tems (T: temperature, fce-Boltzmann const.) [1^ and can be absorbed into a redefinition 
of e, K and 7 by e = e/Fx, ... . A growth trial is performed only if it does not result in too 
strong cell deformations, i.e. if excluded- volume interactions are irrelevant (i.e., if dij > 2R). 
Between two deformation trials a cell performs translation and rotation trials. A cell not 
subject to excluded volume interactions needs 2R/6a growth steps of average size Sa to once 
traverse the cell cycle. Hence its average cycle time is r = {2RAt)/{Sarg). Due to rejected 
growth trials the real average cell cycle time is r^j > r if excluded-volume interactions occur. 
In each case a large Vg corresponds to a large cell growth rate and vice versa. 
In 2(i-space one may distinguish between two idealized, extremal situations, a cell configura- 
tion with fixed ends at xi, X2 (Fig. 2a) and a closed arrangement (Fig. 2b, c). As U > X2 — X1 
{U: arc-length of the cell chain), the stretched configuration becomes curved because the 
minimum free energy configuration as well as the non-equilibrium configurations (Fig. 2a) be- 
come curved. As examples of curved configurations we focus in the following on closed and 
hence intrinsically curved cell geometries (Fig. 2b, c). The spontaneous curvature Cq can be 
dropped in Irf-closed geometries ([^ and below). In each simulation we start with a circular 
configuration of eight cells. For small times t the circumference U and the total number of 
cells grow exponentially fast U cc N cc exp(At) with A = (rZn(2))~^ (Fig.3a) [jl^. Due 
to the excluded volume effect exponential cell division is accompanied by exponentially fast 
depletion of free space on time scales t > r. Hence to maintain exponential growth at still cir- 



4 




FIG. 3: Circumference U in units of 2R + S/2 vers, t/r on lin-log (a) and log-log-scale (b). (a) 
the initial exponential regime is the longer the smaller are k and (squares: oc 2*/'^). For small 
K and Vg a crossover occurs in the geometry (e.g. A:r~^ = 500, k = 20, B: k = 100 in (a)), for 
large k first in the growth law changing from an exponential to a power law (oc with a ~ 0.5; 
tr > t) (b)- For sufficiently large times t, U is bent upwards indicating a second crossover to a 
non-circular configuration. The domain size at the crossover is the smaller, the smaller are and 
K (cf. arrows in (b)). 

cular geometry the migration of each cell i into the direction of the locally outward pointing 
normal nj Vi = 1, .., would have to be exponentially fast, too. The driving force Fj for the 
migration is mainly given by the repulsive part of V^j'^ since due to the elongation of the cell 
axis and the introduction of new cells, the cell-cell-distance is reduced to dij < 2R + 6/2. Its 
component into the direction of is F„- = Fn oc P oc 1/r foi r R thus decreases with 
r (for circular configurations we can drop the index i here) . Since the migration dynamics is 
strongly dissipative, the force does not suffice to relax the cell-cell distance to its equilibrium 
value between two successive cell divisions. For small k, the configuration then roughens at 
still exponential growth (scenario (1), Fig. 3a). For large k the cell-cell distance shrinks until 
it reaches the lower cutoff-value dij = 2R. After this happens only those migration trials 
which fall into the angle interval /3 oc 1/r are successful, others would result in a forbidden 
large deformation of the cells (i.e., dij < 2R). From c/r/rft oc F„ oc 1/r we expect a change 
of the growth law to U oc r oc ^/t at still perfectly circular cell configurations. This is indeed 
found within some limits from the crossover effects (scenario (2), Fig. 3b). For large times 
a second crossover occurs at which the cell configuration starts to deviate significantly from 
a perfect circle, indicated as a characteristic bend upwards in the growth curve (Fig. 3b). 
The instabilities occur at typical domain sizes that grow with (Fig. 2a; arrows in Fig. 
3b) and k, (Fig. 2c, points (A) and (B) in Fig. 3a). 
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The nature of the geometric instabihty can best be understood in a simphfied continuum 
modeL Connecting the middle-points of the cells in a given configuration we arrive at a 
closed curve in 2(i-space. As in ref. we parameterize the curve by the position vector 
r{a,t) where a G [0,1) is a parameter, t the time. r{a,t) and its first two derivatives are 
periodic in a. Note that U = Jq ^/gda where g = dar{a, t)dar{a, t) is the determinant of the 
metric tensor. The mechanisms that contribute to the dynamics are (a) local proliferation of 
arc-length due to cell divisions, (b) stabilization of stretched structures by a bending energy, 
(c) constant (at least on the average) circumference if cell division is switched off. We only 
consider cell configurations near the instability where intersections of the curve with itself 
cannot occur. Then, 



1 dr{a,t) 



/i dt 



(1) 



fi is the mobility. The prefactor ensures reparametrization invariance. JF = J^Q[a,t) — 

Jq A{a, t)y/gda where A(a, t) is a Lagrangian multiplier field which ensures that the condition 
for the local proliferation of length, that is specified below, is fulfilled. JFq = k/2 /Q^(c(a, t) — 
coy^/gda is the bending energy (c(a,t) is the local, cq the spontaneous curvature). To get 
a closed set of equations, we now have to fix a condition for the local metric by ^t^/g = 
daZdadtZ/y/g = f- f dcscribcs the growth of the metric and has to be given externally. 
With the Serret-Frenet relations |]l8l we find fiy/g{cWn + dgWt + c^A - d^A) = f with Wn = 
-{l/^){SJ^o/Sr)n = {K/2){c^ + 2d',c-ccl) {if d^co = 0) and Wt = -{1/ ^){6J^o/6r)t = 0. 
f has been chosen to / = fEGR = '^\f9 in the exponential growth regime (EGR) and 
/ = fpGR = o'y/gc^ in the power-law growth regime (PGR), a measures the strength of 
growth. The choice of fEGR insures exponential growth of the metric, the choice of fpcR 
can be shown to correspond to a normal growth velocity oc c (oc /5 for r ^ R, cf. Fig. 
lb) and to a transversal velocity Gt that obeys the relation dsGt = (from a transformation 
A = A + Gn/c and Gt = —dsiGn/c)). Note that cq can be absorbed into a redefinition of 
A = A — (/t/2)co. These settings result in the homogeneous eqns. dtr = or =^ r oc e'^* 
for EGR and dtr = cr/r ^ r oc \/t for PGR in accordance with the computer simulations. 
Performing a linear stability analysis around a homogeneous growing circle (i.e., in polar 
coordinates: r(a, t) = ro(t) + exp(w(g)t + i2nqa), A{a, t) = Ao(t) -|- rjQ exp{w{q)t + i2'7iqa) 
where r determines the dynamics of A), we find 



+ for EGR 



r4(l+g2) 



^ow^ /or PGR 



'0 



w{q = 0) = cr > for EGR, so homogeneous deviations grow. For PGR w{0) = —a/r^, so 
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homogeneous perturbations are damped out. w{q = 1) = which imphes that a translation 



of the circle is marginally stable ||T6|. A second zero is at = yX(2) for EGR and Qc = 
1+X{1) + Jl + 6X(1) +X(l)2)/2 for PGR where X{z) = r^'a/ifin). For g ^ oo, 



w{q) < i.e., short wavelength perturbations are damped out. For a = the circle is the 



only stable solution. The fastest growing mode is ~ ^X{z)/2 for ^ 1- and qm 
grow with increasing X{z), i.e. the domain sizes at the geometrical instabilities decrease with 
increasing growth rate and decreasing bending rigidity. These results are in agreement with 
the tendencies found in the computer simulations. From the analytical treatment it becomes 
clear that the instability occurs at domain sizes where the bending energy does not suffice to 
smooth the roughening effect of the cell proliferation. An increase of the proliferation rate 
has the same effect: the "time" between two subsequent cell division doesn't suffice anymore 
to smoothen local undulations. For PGR, cells at positions with larger than average local 
curvature have a smaller than average mitotic cycle thus divide faster than cells at positions 
with smaller local curvatures. This further increases the difference in the local curvatures 
in a self-enhancing process resulting in a folding of the domain. For EGR, the increment 
in the metric is proportional to the metric already proliferated. The folded structures are 
not in equilibrium: for the l(i-ring structures, the undulations regress if the cell division is 
suppressed. There is no contribution due to shear. 

Also 2(i-buds formed in 3d may reorganize into a perfectly smooth layer on large time scales 
and reduce shear energy. This line of argument is supported by observations which show 
that cell assemblies can behave as viscoelastic fluids: under compression between two plates, 
cell assemblies first deform followed by a reorganization of the cells in order to reduce elastic 



energy contributions |jT9|. A separate experiment identified cells to diffuse even in aggregates 
11^ . However, growth in a 2(i-surface for some cases is expected to be at least temporarily 



accompanied by shear stress that elevates the energy barrier for the instability, but is not 
expected to change its nature. 

The folding principle neither depends on details of the cell division algorithm nor on the 



particular shape of the NN-potential suggesting some universal feature Hence this 



mechanism may be present in all one layered epithelial tissues including those embedded 
in soft connective tissue during development and the maintenance of tissue if the rate of 
proliferation is sufficiently large either by internal or external stimuli. A typical biological 
example for the folding mechanism may be the post-irradational situation in crypts 



Crypts form pear-shaped pockets in the intestinal wall [^. They are responsible for the 
maintenance of the high cell turnover in the intestine. After irradiation the cell proliferation 
increases and a fingering-like instability occurs in their fioor or walls. Folding is also observed 
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in different geometries such as the skin after treatment with growth factors or during the 



molting and growth of the epidermis of arthropods under hormonal control p2 |. 

For a particular biological situation all simulation parameters can be related to experimental 

quantities p3| . 

DD. is very grateful to M. LoefHer for stimulating discussions and bringing refs. 
into his attention. He further thanks G. Forgacs, M. Kschischo, CS. Potten, U. Seifert and 
especially G. Gompper for discussions or suggestions to the manuscript. 
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